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Abstract. We describe the first discrete-time 4-dimensional numerical application 
of Regge calculus. The spacetime is represented as a complex of 4-dimensional 
simplices, and the geometry interior to each 4-simplex is flat Minkowski spacetime. 
This simplicial spacetime is constructed so as to be foliated with a one parameter 
family of spacelike hypersurfaces built of tetrahedra. We implement a novel two-surface 
initial-data prescription for Regge calculus, and provide the first fully 4-dimensional 
application of an implicit decoupled evolution scheme (the “Sorkin evolution scheme”). 
We benchmark this code on the Kasner cosmology — a cosmology which embodies 
generic features of the collapse of many cosmological models. We (1) reproduce 
the continuum solution with a fractional error in the 3-volume of 10“® after 10000 
evolution steps, (2) demonstrate stable evolution, (3) preserve the standard deviation 
of spatial homogeneity to less than 10“^° and (4) explicitly display the existence 
of diffeomorphism freedom in Regge calculus. We also present the second-order 
convergence properties of the solution to the continuum. 


PACS numbers: 04.20.-q, 04.25.Dm, 04.60.Nc. 

1. Regge calculus as an independent tool in general relativity 

In this paper we describe the hrst fully (3 + l)-dimensional application of Regge calculus 
i0 to general relativity. We develop an initial-value prescription based on the standard 
York formalism, and implement a 4-stage parallel evolution algorithm. We benchmark 
these on the Kasner cosmological model. 

We present three hndings. First, that the Regge solution exhibits second-order 
convergence of the physical variables to the continuum Kasner solution. Secondly, 
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Regge calculus appears to have a complete diffeomorphic structure, in that we are 
free to specify three shift and one lapse condition per vertex. Furthermore, the four 
corresponding constraint equations are conserved, to within a controllable tolerance, 
throughout the evolution. Finally, the recently-developed decoupled parallel evolution 
scheme (the “Sorkin evolution scheme”) yields stable evolution. 

Although we have taken just the first few steps in developing a numerical Regge 
calculus programme, every indication (both numerical and analytic) suggests that it 
will be a valuable tool in the study of gravity. Our numerical studies, together with 
analytic results [|] should put to rest some of the recent concerns about Regge calculus 
i, 1,0 — it does appear to be a viable approximation to general relativity. 

Einstein described gravitation through the curvature of a pseudo-Riemannian 
manifold. Regge calculus, on the other hand, describes gravity through the curvature 
of a piecewise flat simplicial pseudo-Riemannian manifold. The fundamental platform 
for Regge calculus is a lattice spacetime, wherein each lattice cell is a simplex endowed 
with a flat Minkowski geometry p|, ||, ||. The physical and geometric basis of Regge 
calculus distinguishes it from all other discretizations of general relativity. One applies 
the principles of Einstein’s theory directly to the simplicial geometry in order to form the 
curvature, action and field equations . This is in stark contrast to the finite difference 
approach, where one starts with a representation of the continuum held equations and 
proceeds to discretize them over a grid of points. The goal of our work is to evaluate 
the relative strengths and weaknesses of Regge calculus. 

The geometrically transparent nature of Regge calculus should be useful in the 
interpretation of simulations. One example of this is the Kirchoff-like form of the 
contracted Bianchi identity BO in Regge calculus. All edges in the lattice geometry, 
in a strict sense, carry a how of energy-momentum, and the sum of these hows at each 
vertex is zero. This single example illustrates that every term in every Regge equation 
has a clear geometric interpretation. 

We choose to benchmark our code on the Kasner cosmology because it has a well 
dehned solution with symmetries. More importantly, however, the Kasner model is 
a prototype of the Belinsky-Khalatnikov-Lifschitz mixmaster oscillation generic to all 
crunch cosmologies []^|. Unlike all previous Regge simulations, we impose no symmetries 
on the model — it is free to exhibit all dynamical degrees of freedom. The symmetries 
of the Kasner cosmology are encoded in the freely specihable portion of the initial data, 
and the full Regge initial-value problem is then solved. The solution obtained displays 
encouraging agreement with the continuum Kasner solution, and this agreement is found 
to be remarkably well preserved during evolution. The code used in these calculations 
is applicable to x spacetimes, as well as the x topology presented here. 

There are two major areas of active research in numerical relativity today. The 
Erst effort deals with gravitational wave generation by astrophysical processes. The 
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Laser Interferometric Gravity Wave Observatory, which should be operational within a 
few years, will require a set of gravity wave templates from colliding and coalescing 
black holes. The second area of research deals with the structure of cosmological 
singularities in inhomogeneous cosmologies. The gravity wave problem is dominated 
in many ways by boundary conditions, whilst the cosmological singularity problems 
have no boundary. We have chosen to focus our initial application of Regge calculus 
toward issues in inhomogeneous cosmology. 

In this paper we will develop the simplicial lattice used in the calculations (section 
H), and then briefly describe the geometric and dynamical structure of Regge calculus 
(section These ideas will be applied to motivate a York-style two-surface initial-value 
prescription for Regge calculus (section H), which is then specialized to the Kasner 
cosmology (section §). With initial data in hand, we proceed to the time evolution 
problem in section P, and apply a four-step implementation of the Sorkin evolution 
scheme. We examine the convergence of the Regge solution to the continuum, and 
conhrm the existence of simplicial gauge freedoms by demonstrating that the constraint 
equations are satisfactorily preserved during evolution. Finally, in section |^, we discuss 
our future plans for numerical Regge calculus. 

2. Kinematics and the Quantity Production Lattice 

The hrst issue one must address when beginning a Regge calculus simulation is the 
choice of lattice structure. How do we choose to represent a spacetime geometry with a 
lattice geometry? We are aided in our decision by the following four guiding principles 

(i) We want the simplicial spacetime foliated into 3-dimensional tetrahedral spacelike 
hypersurfaces. Although these hypersurfaces will be geometrically distinct, their 
lattice structure should be identical. The spacetime sandwiched between any two 
of these surfaces should be decomposed into simplexes. 

(ii) Each spacelike hypersurface should have the topology of a 3-torus (T^). 

(iii) The simplicial structure sandwiched between two adjacent hypersurfaces must be 
consistent with the recently proposed Sorkin evolution scheme 0. 

(iv) The lattice connectivity, or local topology, of the 3-geometries should be maximally 
homogeneous. In other words, the lattice structure at one vertex should be identical 
to all others. 

These four guiding principles have led us to the choice of lattice described below. It is 
in this sense that we have introduced a (3 -|- 1) split of spacetime in Regge calculus. 

The fundamental geometric element in Regge calculus is the four-simplex, consisting 
of hve vertices, ten edges, ten triangular hinges, and hve tetrahedra. In this section we 
build the four dimensional simplicial lattice by constructing a simplicial three geometry, 
and then carry the triangulation of forward in time, to obtain a simplicial manifold 
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with topology x 

The best tetrahedral subdivision of that we know is based on isosceles tetrahedra, 
and will be referred to as the Quantity Production Lattice (QPL) |^. A major 
advantage of this lattice is that it is easily refined whilst maintaining local homogeneity, 
in the sense that the connectivity at each vertex is the same. The QPL consists of 
isosceles tetrahedra only in flat Euclidean space, with slight variations in the edge lengths 
introducing curvature about edges in the three-lattice. 

The QPL may be constructed from a single cube. Begin by identifying opposing 
faces, thus fixing the global topology. Subdivide the cube into smaller cubes, until the 
desired resolution is obtained, and then introduce a new vertex at the centre of each 
small cube. Joining these centred vertices together yields a new cubic lattice which 
pierces the faces of the original. Finally, join the centred vertices to each of the eight 
vertices of the original cube in which they reside, creating four tetrahedra through each 
face of the original cubic lattice. This completes the construction of the QPL. Figure |I] 
displays some of the local structure of the resulting lattice. 



B A 



Figure 1. Two different views of the Quantity Production Lattice, a) The 
vertex data structure, which when carried from vertex to vertex generates the 
three-geometry. Several A-, B-, C-, and D-type vertices are also shown, and 
these are defined in the text below, b) The local lattice structure about a single 
leg. The four tetrahedra piercing a single face of the original cubic lattice are 
shown “exploded off”. 
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Fourteen edges emanate from each and every vertex in the QPL. Six legs lie along 
cube edges, and 8 diagonal edges join the two cubic lattices together. These two types 
of edges will be refered to as cube-aligned and diagonal edges, respectively. The QPL, 
as outlined here, is isomorphic to the right-tetrahedral lattice, which can be formed 
by adding face and body diagonals to a cubic grid. Despite this, we have found that 
assigning edge lengths using the right tetrahedral approach leads to minor numerical 
difficulties, due to degeneracy in the hrst derivatives of the Regge equations. This has 
perhaps impeded previous work with (3 -|- 1) formulations of Regge calculus M. 

We now have in hand the simplicial three geometry, and turn to the construction 
of the four geometry in which this initial slice is to be embedded. The four geometry 
is constructed by dragging vertices forward, one by one, until the whole initial surface 
has been replicated, and the intervening region hlled with four-simplices. We refer to 


the resulting lattice as a Sorkin triangulation |15, 16 


Another useful approach to the construction of the four geometry is the “vertex data 
structure”, whereby each vertex in the lattice has, in so far as is possible, an identical 
local structure. In general, the Sorkin approach to the evolution problem is inconsistent 
with such a lattice. The major advantage of Sorkin triangulations over a four geometry 
built using a vertex data structure is that it allows the use of an order Ny (number of 
vertices in each three-geometry), parallel, four-step evolution algorithm. That is, one 
quarter of the vertices in the three-geometry may be evolved together, in parallel. 

We construct the Sorkin triangulation as follows. Identify four classes of vertices 
in the QPL, two types on the original cubic structure, and the remaining two on the 
cube-centred lattice . The vertex types are dehned as 


A-type: Half the vertices on the original cubic lattice, sharing no common edges. 

B-type: The remaining vertices on the original cubic lattice. 

C-type: Half the vertices on the cube-centred lattice, not sharing common edges. 

D-type: Remaining vertices on the cube-centred lattice. 

These identihcations are indicated in hgure |I]. We now construct the four geometry 
from the simplicial three geometry. 

Begin by selecting a single A-type vertex, say, and drag it forward in time, 
to create a new vertex A^. In the process we introduce 14 brace edges over the 14 
legs emanating from Ak in the spacelike hypersurface (H;), together with a vertical 
(timelike) leg joining A^ to A^. This timelike edge is taken small enough to ensure that 
the brace edges are spacelike. Each brace edge joins the newly evolved vertex A^ to 
the corresponding vertices in the previous slice which share an edge with A^. Figure ^ 
shows the situation in (1 -|- 1) dimensions. From the hgure it is clear that “evolving” 
an A-type vertex creates a 2-simplex (triangle) over each 1-simplex (edge) meeting at 
Afc. In (3 -|- 1) dimensions, the procedure creates a 4-simplex above each tetrahedron 
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Figure 2. The {1 + 1)-dimensional analogue of the Sorkin scheme, applied to a 
region of the lattice, a) On the initial hypersurface, we identify A- and B-type 
vertices, the only requirement being that no two vertices of a given type share an 
edge. In (3 + 1) dimensions, the lattice breaks into four classes of vertices under 
this requirement, b) All A-type vertices are carried forward together, creating 
(n+l)-simplices “above” each n-simplex in the initial slice, c) The B-type 
vertices are dragged forward. This completes the (1 + 1)-+) evolution procedure, 
replicating the original hypersurface and triangulating the intervening region. 


(3-simplex) which shares the vertex Ak. Having identihed the 4-simphces in this way, 
we may tabulate the new tetrahedra and triangles which they contain. This process can 
be repeated on all A-type vertices simultaneously, since none of these vertices share a 
common edge. 

Now consider a single B-type vertex. Dragging it forward yields 8 brace legs between 
slices, one over each of the 8 diagonal legs in The original vertex is also joined to 
6 A-type vertices which have already been dragged forward, so joining the new vertex 
to these creates 6 spatial legs on S^+i. A vertical edge joining the new vertex to its 
counterpart on the original slice is also created. Again, all B-type vertices can be 
dragged forward in parallel, and this process creates a set of 4-simphces in a manner 
similar to the A-type evolution step. 

Next the cube-centred C-type vertices. Each of these creates 6 brace legs between 
slices, since each C-type vertex in is connected to 6 D-type vertices, which have yet to 
be carried forward. Also created are 8 spatial diagonal legs on S;+i, joining the C-type 
to surrounding A- and B-type vertices, together with a single timelike edge. Finally, 
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the D-type vertices are dragged forward, creating 14 new spatial legs on together 
with a single timelike leg. 

In this way the initial snrface is carried forward in time, constructing an identical 
three geometry on S;+i, and hlling the intervening region with four-simplices. The 
resulting (3 + 1) Regge spacetime is consistent with our four guiding principles. 

3. Dynamical structure of Regge calculus 

In section ^ we developed the kinematic structure of the simplicial spacetime. We now 
apply the Regge form of Einstein’s held equations to the lattice. 

Regge [0,1^ constructed the simplicial form of the Hilbert action, and obtained what 
have become known as the Regge equations via a variational principle. The independent 
variables in this approach are the lattice legs. Each edge in the lattice has associated 
with it a simplicial equivalent of the Einstein equations, which take the form Q 

2L 9Ah _ 

T/* 2 -^ 

h 

where the summation is over all triangles h which hinge on the leg Ij. The defect eh 
about h is the angle that a vector rotates when parallel transported around the hinge. 
Ah is the area of the hinge, and VI is the 3-volume of the Voronoi cell dual to edge L. 
Equation ^ is the Regge calculus analogue of = 0. 

When applied to our lattice, there are three types of Regge equation, depending 
on the type of edge in the lattice it is associated with. The equations associated with 
edges in a spacelike hypersurface may be viewed as “evolution”-type equations, since 
they couple information on that surface to the edges lying above and below. This is the 
Regge equivalent of the second order evolution equations in the ADM (3 -|- 1) split [^j. 

Equations associated with brace edges, those lying between successive spacelike 
hypersurfaces, only involve quantities on and between the two slices. In this sense they 
are “constraint” or first order equations. Similarly, the “timelike Regge equations” 
obtained by varying the vertical edges between two spacelike hypersurfaces are first 
order equations. 

The Einstein tensor in the continuum is a rank 2 symmetric tensor, and hence has 
two indices. E. Cartan has provided us with a geometric interpretation of this tensor 
and its components. He showed that the Einstein tensor is expressible as a double dual 
of the sum of moments of rotation [Q over the 2-dimensional faces of a 3-volume. Thus 
one index can be interpreted as the orientation of that 3-volume, and the other can be 
associated with the orientation of the dual of the sum of moments of rotation. 

The Cartan approach provides a derivation of the Einstein equations independent 
of the action, and has been successfully applied to Regge calculus [Q. The usual 
Regge equations are recovered as a sum of moments of rotations. Thus for each edge 



L in the Regge spacetime the voronoi 3-volume V£ dual to the edge, and each of 
the moment of rotation vectors over the 2-dimensional faces of are parallel and 
directed along L, yielding a single Regge equation per edge. It is in this sense that the 
Regge equation is the double projection of the Einstein tensor along edge L (ie: Gll)- 
Since the timelike edges in our Kasner simulation are orthogonal to the homogeneous 
spacelike hypersurfaces, the Regge equations associated with these edges (the “timelike 
Regge equations”) are the canonical Hamiltonian constraint equations. Similarly, the 
other hrst-order Regge equations, associated with the brace edges, are the only Regge 
equations that carry components of the momentum constraints (Goi)- Thus we use three 
of these brace equations per vertex as momentum constraints, and the timelike Regge 
equation as the Hamiltonian constraint. 

In continuum general relativity, the ten Einstein equations are functions of the ten 
components of the metric tensor. However, of these ten metric components, only six 
functions per spacetime point may be considered truly independent, since we have the 
freedom to choose a system of co-ordinates. Correspondingly, there are four relations 
per point amongst the Einstein equations themselves, the contracted Bianchi identities 


A similar structure is found in Regge geometrodynamics. Each vertex in the 
lattice may be associated with a unique set of N edges, and each of these edges has a 
corresponding Regge equation. However, the simplicial form of the contracted Bianchi 
identity holds at each vertex |]^ , a total of four constraints per vertex. So there can only 
be A^ —4 truly independent edges per vertex, with the remaining four edges representing 
simplicial gauge freedom. This corresponds to the lapse and shift freedom in continuum 
(3 -I- 1) relativity |Q]. The four redundant Regge equations per vertex introduced in 
this process may be compared with the four constraint equations in continuum (3 -|- 1) 
general relativity. 

The hnite rotations involved in the simplicial contracted Bianchi identity do not 
commute, so these “identities” are only approximate 0. In the inhnitesimal limit, 
however, it is expected that the Regge Bianchi identities become exact in correspondence 
to the continuum. This failure to precisely conserve the constraints is no worse than the 
failure of a hnite difference scheme to exactly conserve the numerical constraints. We 
demonstrate this conservation of energy-momentum in section p. 


4. A York-type initial-value prescription for Regge calculus 

Before evolution can proceed, we must construct initial data consistent with the Regge 
constraint equations. All previous initial data constructed in Regge calculus has been 
for the special case of a moment of time symmetry, where the problem reduces to the 
requirement that the scalar curvature of the initial three geometry is zero |^. Successful 
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time-symmetric Regge calculations include initial data for single and multiple black holes 
1^, 115 , Friedmann-Robertson-Walker and Taub cosmologies ^,1^, and Brill waves 


on both flat and black hole backgrounds. 

The major simplihcation of calculating initial data at a moment of time symmetry 
is that the problem reduces to the calculation of the purely three-dimensional scalar 
curvature, and does not involve the four dimensional lattice structure |^. In 
spacetimes of astrophysical interest there will not ordinarily be a moment of time 
symmetry to simplify the initial-value problem. Indeed, this is the case in the Kasner 
class of cosmologies considered here. The discrete time Regge lattice neccesitates 
the construction of fully four-dimensional, two-slice initial data (a “thin sandwich 
approach”). We apply a novel two-slice initial-value formalism, to be described in 
more detail elsewhere il , to construct initial data for the Kasner cosmology in Regge 
calculus. 

The two-surface approach outlined below is constructed to mirror, in as far as 
possible, the York ^ technique for constructing initial data in the continuum, rather 
than the Belasco-Ohanian ||^ two-surface formulation. The advantage of the York 
approach is that it clearly delineates the true degrees of freedom of the gravitational 
held from the embedding quantities. York begins by conformally decomposing the three- 
metric, 

introducing the conformal factor ifj, and base metric We are then free to specify 
the 


• base metric (“Where”). 

• momenta (“How Fast”), and 

• The trace of the extrinsic curvature, TtK (“When”). 

Once these have been hxed, the only remaining variables are the conformal factor 
and the gravitomagnetic three-vector potential W\ The four constraint equations are 
used to calculate these quantities. 

The key to our simplicial two-surface formulation of York’s procedure is the 
simplicial representation of the dynamical degrees of freedom and embedding variables. 
The identihcation presented below is by no means unique. In other words, we provide 
a representation of the degrees of freedom and not a diffeomorphically invariant, 
conformally invariant representation of the dynamic degrees of freedom. To this end, it 
is convenient to perform a conformal decomposition on each hypersurface, yielding leg 
lengths in (k = 0,1) of the form 

hj = 'i’ij hj 

where hj is the base leg between vertices i and j lying in one hypersurface, the simplicial 
equivalent of the base metric The conformal factor is defined on the vertices of 
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the 3-lattice S, and is applied to the edge between vertices i and j using a centred, 
second order, approximation 

= \i.A + i’j) ■ 

Within each hypersurface, 14 edges emanate from each vertex; 6 cube-aligned, 
together with 8 diagonal legs. It is convenient at this stage to view the lattice, and 
the three-lattice in particular, from the viewpoint of a vertex data structure. The three 
geometry can be constructed by applying the generator shown in figure to each vertex, 
which consists of 3 cube-aligned legs, together with 4 diagonal edges. The four geometry 
sandwiched between the two hypersurfaces contains a single brace leg over each of these 
7 edges per vertex, and so a similar decomposition into edges per vertex can be obtained 
for the brace legs. However, the exact structure of these braces depends on the choice 
made in the construction of the four lattice, and in particular, each type of vertex 
(A,B,C,D) will have a different arrangement of braces about it. 

Since we are dealing with a full four-dimensional region of the lattice, the initial- 
value problem requires that we specify lapse and shift, otherwise we would be unable to 
construct the thin-sandwich analogue of TtK and These quantities both depend 
on derivatives of the 3-metric, the lapse function and the shift vector. Our freedom 
to freely choose lapse and shift is linked to the simplicial contracted Bianchi identities, 
discussed above [Q . Although Regge calculus deals directly with geometric, co-ordinate 
independent quantities, we must still specify how vertices on the initial three geometry 
are pushed forward in time. 

The York prescription outlined above allows us to freely specify the base three- 
geometry (the “where” |^), TrA' (“when”), and the momentum conjugate to the true 
dynamical degrees of freedom of the gravitational held (“how fast”). 

A two-surface formulation of these ideas requires that the momenta terms be split 
across both surfaces. In direct analogy to the assignment of the base three-metric, we 
specify freely all base legs on the initial hypersurface Eq. To avoid specifying everything 
on either surface (and thus obtaining a Belasco-Ohanian style initial data prescription 
@1), we must not specify the conformal factor -0 on Sq. The conformal factor at each 
vertex on Ei is freely specihed instead, in analogy to the hxing of TtK. The hnal step 
is to identify the Regge analogue of We choose the 4 diagonal base legs per vertex 
on El to represent this freedom, since, together with the lapse and shift choice, they 
enable us to specify part of the change in the base “metric” across the sandwich. 

Our two-surface formulation on the lattice, in the spirit of the standard York 
decomposition, may be summarized as follows. Specify 


• All 7 base edges kj per vertex on Eq (“Where”), 

• The 4 diagonal base edges per vertex on Ei (“How Fast”), 

• The conformal factor ipij at each vertex on Ei (“When”), and 
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• 4 lapse and shift conditions per vertex (“How Fast” and “When”), 

which leaves the conformal factor -0 on Eq, the 3 spatial cnbe-aligned edges per vertex on 
El, and the 4 braces per vertex which lie over diagonal edges as the true representations 
of the geometric degrees of freedom. These 8 unknowns per vertex may be calculated 
by solving the 8 Regge constraint equations per vertex which are available during the 
solution of the initial-value problem. The constraints are associated with the 7 braces 
per vertex, together with the timelike edge joining the vertex on Eq to its counterpart 
on El. 

It is possible to introduce sophisticated definitions of lapse and shift by mirroring the 
continuum structure [^, however a more natural approach is to select appropriate legs 
in the lattice, and define them to represent the simplicial gauge freedom. In this spirit, 
we dehne lapse to be the proper time measured along a timelike leg joining a vertex to 
its counterpart on the next hypersurface, and the shift to be a combination of three legs 
that uniquely fix the location of the new vertex above the previous hypersurface. 

Detailed definitions of shift differ as to which type of vertex (A,B,C or D) we are 
considering. The evolution of an A-type vertex creates 14 brace edges stretching between 
E; and the nascent hypersurface E^+i, together with a timelike “Regge lapse” edge. Six 
of these edges will lie above cube-aligned edges in E/, and it is these braces that we select 
as our simplicial shift freedom. In particular, at a given A-type vertex, the shift edges 
correspond to the brace legs above the 3 cube-aligned edges in the vertex generator 
shown in figure ^ This definition of shift is also applied to the C-type vertices. 

In the case of B- and D-type vertices, assigning the lengths of these braces will 
not determine the position of the evolved vertex. This can be seen in figure where 
assigning all brace lengths will only rigidify the “A”-type vertices, leaving the “B”-type 
to flap. This would result in an ill-conditioned system of equations, since the location 
of the B-type vertices is completely undetermined. This (1 -|- l)-dimensional example 
suggests that we define the shift edges on B- and D-type vertices by assigning the lengths 
of the 3 cube-aligned edges in E^+i. The three edges chosen are always those associated 
with the vertex-by-vertex generator shown in figure |I]. It is clear from this discussion 
that care must be taken in choosing the simplicial counterparts of lapse and shift. 


5. Thin-sandwich initial-value approach applied to the Kasner cosmology 


We now apply the initial-value formalism outlined above to the Kasner class of 
cosmologies. The metric for these x vacuum solutions of Einstein’s equations 
takes the form 


ds^ = — 


de + + fP^dy^ + t^P^dz^ 


( 2 ) 
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where the unknown constants Pi satisfy 

Pi+P2+Pz=p\+pI+pI = 1- (3) 

We construct the two surface initial data such that the initial surface, Sq, is at f = 1, 
and El is at f = 1 + At. The base leg lengths on Sq are assigned the flat-space values 
obtained by setting t = 1 in (^. This yields 

= Axk (4) 

for the three cube-aligned base edges per vertex on Eg, and 

for the 4 diagonal base edges per vertex. 

The conformal factor on Ei is the simplicial equivalent of setting Tri^, and we choose 
to take 


-0 = 1 on El. (6) 

The kinematic degrees of freedom correspond to the 4 diagonal legs per vertex on Ei, 
and these are set in the same manner as the base diagonal legs on the initial slice, using 
(I) evaluated at t = 1 -|- At. This yields 

[ 2 . = 1 (t^p^Ax^ + + e^^Az^) . ( 7 ) 

The only edges which remain to be hxed are our choice of shift and lapse. The lapse 
edge is assigned the squared proper time between slices, 

( 8 ) 

and the squared lengths of the shift edges are obtained by applying a power series 
expansion along the continuum geodesics. Accurate to third order in the lattice spacing, 
the expansion between vertices i and j is 

= Qf.uAx^jAx'^j + gc.pTy^AxfjAx^jAx’ij + ... (9) 

where the continuum metric g^u is used to obtain T^^. The spatial edge length 
assignments above are identical to the expressions obtained from (|^), to third order. 
In principle we could ignore the cubic term in the expansion, and take the squared shift 
edge lengths to be 

lij = g^uAx^jAx’^ij, ( 10 ) 

however it was found that the initial-value problem converged faster using the higher 
order expansion. The simplest approximation was tested, and found not to change the 
character of the initial-value solution appreciably. These expansions are used only whilst 
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constructing initial-value data, and only when we wish to compare with a continuum 
metric. 

The shift conditions on A- and C-type vertices are applied to brace legs lying above 
cube-aligned edges, and the geodesic-length approximation (|]) for these yields 

+ PkAxlAt. ( 11 ) 

For the B- and D-type vertices, shift is applied to the three cube-aligned base spatial 
edges per vertex in Si, where the series expansion takes the form 

!l^ = -At^ + (1 + Atf^^ Axl (12) 


In all calculations shown in this paper, we use the QPL obtained by two barycentric 
subdivisions of the original cube. The three geometry consists of 128 vertices, 896 
legs, 1536 triangles and 768 tetrahedra per spatial hypersurface. The region contained 
between two consecutive spacelike hypersurfaces contains 3072 four-simplices. Whilst 
this is a modest grid resolution by modern standards, the major aim of this work is the 
evaluation of Regge calculus as a competitive technique in (3-1-1) numerical relativity; 
we feel that the grid is sufficient for this task, and experiments with increased spatial 
resolution do not effect the conclusions reached below. 

The Kasner exponents pi enter into the freely specihable portion of the initial 
data through equations 0. o and (p!2D. They are not used at any other stage 
of the calculations, and in particular, do not appear in the evolution problem. The 
information borrowed from the continuum solution in this way could, in principle, 
be constructed within the Regge framework by seeking a solution with the required 
symmetries. Analytic calculations in which the Regge lattice is forced to maintain the 
symmetries of the Kasner solution have been undertaken p7| , p8| . These conhrm that 
the Regge equations produce the expected Einstein equations for the Kasner cosmology 
to leading order in the limit of very hne discretizations. Any calculation must begin 
with a suitable ansatz to determine the desired class of solutions; we have chosen the 
freely specihable portion of the initial data to closely model the Kasner cosmology, 
allowing us to compare with the analytic solution. This provides a test-bed for (3-1-1) 
Regge calculus, whilst maintaining maximal freedom in the non-specihable portion of 
the initial data. 

The spatial length scales which appear in the geodesic expansions above must be 
chosen for a particular simulation, noting that the total volume of the initial base 
geometry will be QAAxAyAz. To avoid problems with the Courant condition, the 
timestep is scaled together with the spatial resolution. We choose 


At = min{Ax, Aj/, A^} x Ato, 


(13) 


where Afo < 0.5. The only constraint on the choice of spatial length scale is the 
condition that we require many grid zones across the co-ordinate horizon, which for the 
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metric is located at t, so at t = 1 we require Axj << 1 p9 . 

The initial-value algorithm outlined in section 4 can now be applied to the QPL 
model of the Kasner cosmology. The freely specifiable variables are fixed using equations 
(|)-(|), with equations ([Tl|) and ([T^) being used to enforce the shift conditions. The 
remaining variables at each vertex, namely the conformal factor on Sq, the four diagonal 
brace edges between slices, and the three cube aligned edges in Ei, are calculated by 
solving the available Regge equations. The resulting system of eight equations per vertex 
(a total of 1024 variables for the QPL considered here) is solved using Newton-Raphson 
iteration. 

Finally, we present a sample initial data solution for an axisymmetric Kasner 
cosmology (pi = P 2 = 2/3,pa = —1/3), where we take Ax = Ay = Az = 0.005 
and At = 0.001. The initial-value solution is found to be homogeneous to a very high 
degree. Solving the constraints for the independent variables in this case yields the 
mean values 


= 1.0000000001 
[2 = 0.0000250333 
= 0.0000250333 
!l = 0.0000249833 
ci2 = 0.0000177563 


(14) 


where lx, ly and I 


are the three cube-aligned legs in Si, d is a diagonal brace edge 
between slices, and V’ is the conformal factor on Sq. The standard deviation of both 
the edges and the conformal factor from these mean values is of the order 10“^°, which 
is the convergence tolerance used for the Newton-Raphson iteration. 


6. Time evolution of the Kasner cosmology 


In the previous sections we described the structure of the quantity production lattice, 
and built two-surface initial data for the Kasner cosmology. We now consider the 
evolution of this data. 

The initial data is evolved using the parallel Sorkin evolution scheme described by 
Barrett ef a/ 0. In section ^ we constructed the simplicial lattice from the initial 3- 
geometry by dragging forward individual vertices. This procedure was used because it 
results in a lattice suited to Sorkin evolution. The key realization in the Sorkin approach 
n, Dased on the original examples |^, is that when a vertex is carried forward to 


the next slice, a new Regge equation becomes available below every new edge created. 
Thus a purely local evolution algorithm is possible. The evolution of a single vertex 
in the quantity production lattice creates 15 new edges, and makes available 15 Regge 
equations. 
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In the case of the A-type vertices, 14 brace edges are created between S/ and 
together with the tinielike edge. This completely encloses the region snrrounding the 

14 spatial edges emanating from the vertex in and the timelike edge joining the 
vertex to its fnture connterpart. The Regge eqnations associated with the legs in S/ 
emanating from the vertex A, together with the timelike edge, may be nsed to solve for 
the lengths of the 14 new brace edges and single new timelike edge. A similar sitnation 
applies for the B-type vertices, except here there are 6 new cube-aligned edges in 

8 new spacelike diagonal braces and a timelike leg between the Cauchy slices. In this 
case, the new Regge equations which become available are associated with the timelike 
edge, the 6 cube-aligned braces below the new spatial edge, and the 8 spatial diagonal 
edges in S;. Evolving a C-type vertex creates 6 new brace edges, together with 8 spatial 
edges and a single timelike edge, and the available Regge equations correspond to the 6 
spatial edges in Sq, 8 brace edges, and the timelike edge. Finally, D-type vertices create 
a vertical edge and 14 new spatial edges on the newly formed spacelike hypersurface. 
The evolution equations in this case are associated with the brace edges between slices 
and the timelike edge at the vertex. 

However, this is not the entire story. As we saw in section the simplicial form of 
the contracted Bianchi identities imply that four equations per vertex are dependent on 
the remaining equations. The quantity production lattice has a total of 15 dynamical 
edges per vertex during evolution, but we can only consider 11 of those as independent, 
with the remaining four freely specifiable - the simplicial equivalents of shift and lapse 
freedom. 

We choose to apply the lapse and shift conditions on the same combination of edges 
as in the initial-value problem, however the technique is slightly different. Consider an 
A-type vertex. About this vertex there are 6 brace edges, lying above the cube-aligned 
edges in the current hypersurface. The zero shift condition is applied by demanding that 
opposing brace edges, lying along the same “axis”, have equal lengths. This condition 
is applied to each of the three pairs of such edges. Provided that the homogeneity of 
the initial slice is maintained, this ensures that the timelike edge is locally orthogonal to 
the current hypersurface. Similar conditions are used for the remaining vertices, applied 
either to brace edges (C-type) or spatial edges on the next hypersurface (B- and D-type 
vertices). The lapse freedom is utilised to specify the squared length of the timelike 
edge at each vertex. 

The application of four gauge conditions per vertex leaves 11 dynamical edges and 

15 equations at each vertex in the lattice. Clearly 4 Regge equations per vertex are 
redundant, becoming constraints which may be tracked during evolution. The redundant 
equation associated with lapse freedom is our “Hamiltonian constraint”, the timelike 
Regge equation. When applying shift conditions, we choose the redundant equation 
to lie below the edges upon which the gauge condition is applied. For A- and C- 
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type vertices, shift is applied to brace edges, so the 3 redundant equations per vertex 
correspond to the cube-aligned spatial edges in the vertex generator (see figure |^). For B- 
and D-type vertices, shift is applied to spatial edges in so the redundant equations 
are associated with the brace edges lying below the 3 cube-aligned vertex generator legs 
on S;+i. 

The Sorkin evolution scheme is, by construction, a partially constrained algorithm, 
since the redundant equations associated with shift edges at A- and C-type vertices are 
second order evolution-type equations, whereas at B- and D-type vertices we discard 
first order constraint-type equations. For simplicity we shall refer to the redundant shift 
equations at each vertex as the “Momentum constraints”, although they are in fact a 
mixture of first-order and second-order equations. 

The evolution of an axisymmetric Kasner cosmology (with Kasner exponents 
Pi = P2 = 2/3, ps = —1/3) is shown in figure |], and figure § shows the evolution of 
an expanding fiat space solution, corresponding to the Kasner exponents Pi = P 2 = 
0,p3 = 1. Both runs cover a ten-fold increase in the three-volume. It can be seen that 
the fractional error in the three-volume remains small throughout the evolution. The 
fractional volume discrepancy is defined as 


AK 

'V 


3y 


(15) 


and the summation is over all tetrahedra on a t = constant surface. 

Also apparent, particularly in the axisymmetric Kasner solution in figure |^, are 
high frequency oscillations (typical wavelengths of a few timesteps), which decrease 
in magnitude as the evolution proceeds. It was pointed out to us by R. Matzner 
that gravitational waves in a Kasner background would evolve with such a behaviour 
13^, ^ . We believe we are seeing a grid resolution limit effect reminiscent of such waves; 
however, in our case they are not fully resolvable. The oscillations are controllable, with 
magnitudes varying as and wavelengths reducing as I, where I is a typical edge length 
on the lattice. The presence of these relatively high frequency, low amplitude waves 
does not appear to cause any instabilities, and it is apparent from figures ^ and § that 
the waves gradually die out as the lattice evolves. 

The evolution of the three different classes of spatial cube-aligned edge are shown in 
figures 1^ and and are found to closely match the expected Kasner evolution. In both 
figures the lx and ly edges he atop one another. The standard deviation of the spatial 
edge Iz, shown in part d of the figures, increases gradually throughout the evolution. It 
reaches a few parts in 10“^° after ten thousand timesteps for the axisymmetric Kasner 
solution, and remains below 10“^° in the expanding fiat spacetime of figure 

We now consider the convergence properties of the simplicial cosmological solution. 
Due to the global topology and homogeneity of the solution, we are able to examine 
convergence by reducing the typical scale length of the lattice, whilst keeping the number 
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Figure 3. The axisymmetric simplicial Kasner solution with exponents 
(2/3, 2/3,—1/3). The seale lengths are taken to he Ax = Ay = Az = 0.005, 
and At = 0.001. a) The fraetional error in the spatial 3-volume of the newly- 
evolved spacelike hypersurface is plotted against the proper time measured along 
timelike edges in the lattice. The error represents the fraetional differenee in 
the simplieial 3-volume from the expected linear growth rate, b) Growth of the 

Hamiltonian ( - ) and Momentum ( - ) constraints with time, e) Mean 

value of the spatial legs G ( - ) and G ( - ) against time. The legs display 

the expected power law behaviour, with the ly leg being identieal to G throughout 
the evolution, d) Standard deviation of the spatial leg Iz from mean. The other 
spatial legs lx, ly, and Id behave similarly. 
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Figure 4. The flat Kasner solution, with exponents (0,0,1). The lattice 
maintains homogeneity to a very high degree, with the standard deviation 
remaining below the limit of numerical accuracy (approximately 1 x The 

spatial and temporal scales are Ax = Ay = 0.125, Az = 0.005, and At = 0.002. 
The same quantities are plotted as in the non-flat axisymmetric case shown in 
figure 


of vertices fixed. This is equivalent to looking at a smaller region of the total manifold. 
For the remainder of this section we introduce the scale parameter 6, defined such that 

Ax = Ay = Az = - and At = 0.05 x 5, (16) 

which ensures that the Courant condition remains satisfied, and the total volume of the 
initial base three-geometry is 5^. 

There are several important convergence issues that must be examined in Regge 
calculus, and our fully simplicial (3-1-1) calculation provides an excellent opportunity 
to do so for a specific application. The first, and most important issue to address 
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is the convergence of simplicial solutions to the corresponding solutions of Einstein’s 
equations. The second issue is the consistency of the Regge equations. That is, as the 
typical lattice length scale is reduced, do the redundant Regge equations also converge? 
This latter question is related to the simplicial Bianchi identities, and the conservation 
of energy momentum. We hnd convergence in all solutions and equations. 
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Figure 5. Convergence estimates for the axisymmetric Kasner model shown in 
figure |^. a) The average value of the fractional volume discrepancy, equation 
& evaluated at t = 8, is plotted against the overall scale parameter 5 (A- 
points) . The convergence in the diagonal spatial legs (A-type points), also 
evaluated att = 8, is shown. Both the error in the volume and the diagonal legs 
scales as the second power of 5. This confirms, for the axisymmetric simplicial 
Kasner cosmology at least, that the solution of the Regge equations is a second 
order accurate approximation to the continuum solution. All other spatial legs 
were found to converge as second order in the scale length, b) The average 
magnitude of the constraint equations, evaluated at t = 8, plotted against 
the lattice scale length 6. The momentum constraints (A-type points) display 
slightly better than second order convergence. The Hamiltonian constraints (A- 
type) scale linearly with 5. 

The error in lattice edge lengths compared to corresponding geodesic segments in 
the continuum was examined, together with the fractional difference of the simplicial 
three volume from linear expansion displayed by the continuum. The fractional volume 
discrepancy, (0. was evaluated at f = 8, although due to the fluctuations apparent in 
hgure the values were averaged over several periods of the oscillation. The results are 
shown in hgure ^ 

Figure ^ also shows that the diagonal spatial edges converge to the continuum 
solution second order in the typical lattice scale length. All other spatial edges were 
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found to converge at the same rate. This indicates that the simplicial Kasner solution 
converges to the continuum solution as the second power of the lattice spacing 6. 

Recent work has cast doubt over the convergence of Regge calculus to the continuum 
P, 1^. The calculations of Brewin applied continuum geodesic lengths directly to 
the lattice, and evaluated the Regge equations on the resulting simplicial spacetime. 
Numerical convergence tests indicated that the Regge equations failed to converge as 
the scale length of the lattice was decreased. We have shown explicitly, for one of 
the cases studied by Brewin, that solutions of the Regge equations converge to the 
corresponding solution of the continuum Einstein equations. 

Figure also displays the convergence of the redundant Regge equations as the 
lattice is refined. The standard deviation of the Momentum constraints, a combination 
of both evolution and true constraint-type equations, is shown to converge as at least 
the second power of 6. The Hamiltonian constraint shows linear convergence. The 
convergence of the equations is consistent with the Regge Bianchi identities, and the 
existence of diffeomorphism freedom in Regge calculus. 

The convergence analysis shows that for the Kasner cosmological model, Regge 
calculus is a second order discretization of Einstein’s theory of gravity. 


7. From BKL to Geons; future directions in Regge calculns 


We have successfully performed the hrst fully (3 -|- l)-dimensional calculation in Regge 
calculus, without the imposition of symmetry conditions. In the process we applied a 
novel 2-surface initial value formalism to the lattice, demonstrated the gauge freedom 
implied by the simplicial form of the contracted Bianchi identities, and showed explicitly 
that the solution of the Regge equations converges to its continuum counterpart. 

The simplicial Kasner solution was found to agree well with the analytic 
solution, and maintains a remarkable degree of homogeneity throughout the evolution. 
Convergence analysis showed that all edge lengths converge to their continuum values 
as the second power of the typical lattice length scale, countering recent doubts over 
the convergence of Regge calculus to general relativity. The Regge constraint equations, 
arising from simplicial gauge freedoms, were also found to converge to zero as the lattice 
was rehned. This demonstrates explicitly the existence of gauge freedom in the lattice, 
the simplicial counterpart of coordinate freedom in continuum general relativity. 

Three issues direct our research in the immediate future. First is an analytic 
formulation of the convergence properties of the Regge equations, and their solutions. 
Secondly, we wish to provide an x benchmark of our code, based on a Taub- 
like cosmology. Finally, we are investigating planar numerical perturbations [^] of 
the Kasner cosmology, in order to analyze gravity wave propagation in our simplicial 
spacetime geometry. 
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In the longer term we intend to address three issnes nsing Regge calculns. The hrst 
is related to the wiring of matter terms to the lattice geometry. We cannot think of a 
better approach than to ntilize the Cartan analysis |Q described in section |^. With a 
formulation of matter in Regge calculus, we can begin to study the generic properties of 
collapse in inhomogeneous cosmologies. At both early and late stages in cosmological 
expansion, the effective gravity wave energy density dominates the contributions from 
matter and radiation helds. For this reason we wish to investigate geon states |] of 
the gravitational held 0- We believe that these goals provide a clear direction for the 
future development of Regge calculus. 
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